home *** CD-ROM | disk | FTP | other *** search
/ MacTech 1 to 12 / MacTech-vol-1-12.toast / Source / MacTech® Magazine / Volume 08 - 1992 / 08.03 Jul 92 / Fast Random Numbers / RandomNumbers.a < prev    next >
Encoding:
Text File  |  1991-10-05  |  5.9 KB  |  218 lines  |  [TEXT/MPS ]
  1.         MACHINE MC68020
  2.         MC68881
  3.  
  4. ;-----------------------------------------------------------
  5. ;  WARNING:  These routines require a 68020 or 68030 CPU,
  6. ;  and a 68881 or 68882 floating-point coprocessor!!!
  7. ;-----------------------------------------------------------
  8. ;  This unit implements a "minimal standard" random number
  9. ;  generator using 32-bit integer arithmetic.  This is not
  10. ;  the "best" random number generator possible, but it is
  11. ;  simple, quick, and produces numbers which are "random"
  12. ;  enough for most purposes.
  13. ;
  14. ;  It is based on the "parametric multiplicative linear con-
  15. ;  gruential algorithm" which was originally proposed by D.
  16. ;  H. Lehmer in 1951.  This algorithm generates a sequence
  17. ;  of integers z[1], z[2], z[3]... by repeatedly carrying 
  18. ;  out the following calculation:
  19. ;
  20. ;               z[n+1] = (a * z[n]) mod m
  21. ;
  22. ;  where "a" and "m" are suitable constants.  The successive
  23. ;  "z"s are stored in the variable "seed" in this unit.  
  24. ;  They all lie in the range [0..m-1], so dividing them by
  25. ;  "m" gives real numbers which lie in the range [0.0..1.0].
  26. ;
  27. ;  In 1969, Lewis, Goodman and Miller proposed the choice of
  28. ;  a = 16807 and m = 2**31 - 1 = 2147483647 as particularly
  29. ;  suited for a machine with a 32-bit word length (in their
  30. ;  case, the IBM System/360).  The main problem is that the
  31. ;  product a * z[n] can be more than 32 bits long, leading
  32. ;  to integer overflow on many systems (including the Mac).
  33. ;
  34. ;  In 1979, G. L. Schrage devised an implementation of the
  35. ;  LG&M method which is guaranteed not to produce integer
  36. ;  overflow on any system with a 32-bit word length.  It is
  37. ;  this implementation which is used here, in procedure
  38. ;  "UpdateSeed".
  39. ;
  40. ;  Source:  Stephen K. Park and Keith W. Miller, "Random
  41. ;  Number Generators:  Good Ones Are Hard to Find", Communi-
  42. ;  cations of the ACM, vol. 31, p. 1192  (October 1988).
  43. ;  
  44. ;  This unit was written in MPW Assembler, v3.2.  If the
  45. ;  calling program is written in Pascal, it must use the
  46. ;  interface unit "RandomNumbers.p".
  47. ;
  48. ;                          Jon Bell
  49. ;             Dept. of Physics & Computer Science
  50. ;                    Presbyterian College
  51. ;                      Clinton SC 29325
  52. ;                   CompuServe: #70441,353
  53. ;                         October 1991                       
  54. ;-----------------------------------------------------------
  55. ;  Numerical constants used in the seed-updating algorithm.
  56.  
  57. a        EQU        16807
  58. m        EQU        2147483647
  59. q        EQU        127773    ; m div a
  60. r        EQU        2836        ; m mod a
  61.  
  62. ;  NOTE:  Other possible values for these constants, which
  63. ;  may give "better" results, are:
  64. ;
  65. ;     a = 48271,  q = 44488,  r = 3399  or
  66. ;     a = 69621,  q = 30845,  r = 23902,
  67. ;
  68. ;  keeping m = 2147483647.
  69. ;-----------------------------------------------------------
  70. ;  The random number seed (a 32-bit integer) is a global
  71. ;  variable whose value is preserved between calls to
  72. ;  the random number functions.
  73.  
  74. Seed        DS.L    1
  75.  
  76. ;-----------------------------------------------------------
  77.  
  78. InitRandomSeed    PROC        EXPORT
  79.  
  80. ;- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  81. ;  Initializes the random number seed to a specified value.
  82. ;
  83. ;  Pascal interface:
  84. ;
  85. ;      procedure InitRandomSeed (newSeed : longint);
  86. ;- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  87.  
  88. newSeed    EQU        8
  89. ParamSize    EQU        4
  90. LocalSize    EQU        0
  91.  
  92.         LINK        A6, #LocalSize
  93.         MOVE.L    newSeed(A6), seed(A5)
  94.         UNLK        A6
  95.         RTD        #ParamSize
  96.  
  97.         DC.B        'INITRANDOMSEED'
  98.  
  99. ;-----------------------------------------------------------
  100.  
  101. RandomSeed        FUNC        EXPORT
  102.  
  103. ;- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  104. ;  Returns the current value of the random number seed.
  105. ;
  106. ;  Pascal interface:
  107. ;
  108. ;      function RandomSeed : longint;
  109. ;- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  110.  
  111. result    EQU        8
  112. LocalSize    EQU        0
  113.  
  114.         LINK        A6, #LocalSize
  115.         MOVE.L    Seed(A5), result(A6)
  116.         UNLK        A6
  117.         RTS
  118.  
  119.         DC.B        'RANDOMSEED'
  120.  
  121. ;-----------------------------------------------------------
  122.  
  123. UpdateSeed        PROC
  124.  
  125. ;- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  126. ;  Applies the Lehmer / Lewis, Goodman, Miller / Schrage
  127. ;  algorithm to update the random number seed.  This
  128. ;  procedure is not to be used outside this unit, therefore
  129. ;  it is not EXPORTed.
  130. ;
  131. ;  It stores the new seed in the global variable "Seed",
  132. ;  and leaves a copy in register D1 for use by the calling
  133. ;  routine.
  134. ;- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  135.  
  136. hi        EQU        D0
  137. lo        EQU        D1
  138.  
  139.         MOVE.L    Seed(A5), hi
  140.         TDIVS.L    #q, lo:hi
  141.         MULS.L    #a, lo
  142.         MULS.L    #r, hi
  143.         SUB.L        hi, lo        
  144.         BGT.S        StoreNewSeed
  145.         ADD.L        #m, lo        ; if lo <= 0
  146.  
  147. StoreNewSeed
  148.  
  149.         MOVE.L    lo, Seed(A5)
  150.         RTS
  151.  
  152.         DC.B        'UPDATESEED'
  153.  
  154. ;-----------------------------------------------------------
  155.  
  156. RandomReal        FUNC        EXPORT
  157.  
  158. ;- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  159. ;  Updates the random number seed and returns a real number
  160. ;  in the range (0,1).
  161. ;
  162. ;  Pascal interface:
  163. ;
  164. ;      function RandomReal : extended;
  165. ;- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  166.  
  167. quotient    EQU        FP0
  168. newSeed    EQU        D1
  169. result    EQU        8
  170. LocalSize    EQU        0
  171.  
  172.         LINK        A6, #LocalSize
  173.         JSR        UpdateSeed
  174.         FMOVE.L    newSeed, quotient
  175.         FDIV.L    #m, quotient
  176.         FMOVE.X    quotient, ([result,A6])
  177.         UNLK        A6
  178.         RTS
  179.  
  180.         DC.B        'RANDOMREAL'
  181.  
  182. ;-----------------------------------------------------------
  183.  
  184. RandomInteger    FUNC        EXPORT
  185.  
  186. ;- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  187. ;  Updates the random number seed and returns a longint
  188. ;  in the range [1,max].
  189. ;
  190. ;  Pascal interface:
  191. ;
  192. ;      function RandomInteger (max : longint) : longint;
  193. ;- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  194.  
  195. dividend    EQU        D1
  196. divisor    EQU        D0
  197. result    EQU        8
  198. max        EQU        12
  199. LocalSize    EQU        0
  200. ParamSize    EQU        4
  201.  
  202.         LINK        A6, #LocalSize
  203.         JSR        UpdateSeed
  204.  
  205.         ; The new seed is now in "dividend."
  206.  
  207.         TDIVS.L    max(A6), divisor:dividend
  208.  
  209.         ; "Divisor" now contains the remainder.
  210.  
  211.         ADDQ.L    #1, divisor
  212.         MOVE.L    divisor, result(A6)
  213.         UNLK        A6
  214.         RTD        #ParamSize
  215.  
  216.         DC.B        'RANDOMINTEGER'
  217.  
  218.         END